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A micromechanical model of collapsing quicksand 
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Abstract The discrete element method constitutes a gen- 
eral class of modeling techniques to simulate the micro- 
scopic behavior (i.e. at the particle scale) of granular/soil 
materials. We present a contact dynamics method, ac- 
counting for the cohesive nature of fine powders and soils. 
A modification of the model adjusted to capture the es- 
sential physical processes underlying the dynamics of gen- 
eration and collapse of loose systems is able to simulate 
" quicksand" behavior of a collapsing soil material, in par- 
ticular of a specific type, which we call "living quicksand" . 
We investigate the penetration behavior of an object for 
varying density of the material. We also investigate the 
dynamics of the penetration process, by measuring the 
relation between the driving force and the resulting veloc- 
ity of the intruder, leading to a "power law" behavior with 
exponent 1/2, i.e. a quadratic velocity dependence of the 
drag force on the intruder. 

Key words. Granular matter, Contact Dynamics Simu- 
lations, Distinct Element Method, Quicksand, Collapsible 
soil, Biomaterial 



Introduction 

Despite the ubiquitous appearance of quicksand in adven- 
ture books and movies, its origin and physico-chemical be- 
havior still represent controversial scientific issues in the 
fields of soil and fluid mechanics [J[5J[3J[MS] • ^ has been 
argued repeatedly [S] that, because the density of sludge is 
typically larger than that of water, a person cannot fully 
submerge, and therefore cannot be really "swallowed" by 
any quicksand. 

The fluidization of a soil due to an increase in ground 
water pressure, which in fact is often responsible for catas- 
trophic failures at construction sites, is called by engineers 
the "quick condition" [TIE] , and has been studied exten- 
sively up to now [9lll0l[TTj. Another source of fluidization 
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can be vibrations cither from an engine [12] or through an 
earthquake [3] . While this "liquefaction" or "cyclic mobil- 
ity" phenomenon [T51[T^l[T5] can essentially happen with 
any soil [16] , it is known that samples taken from natural 
quicksand usually show quite specific, but anomalous rhe- 
ology depending on the peculiarities of the material com- 
position and structure [T7] . Also in dry quicksand [TSIH^] a 
fragile/metastable structure leads to interesting material 
behavior like collapse and jet formation after impact of 
an intruder, although the microscopic material properties 
are obviously quite different. Nevertheless, wet and dry 
quicksand are more similar than one would think. Based 
on results of real and in situ quicksand measurements we 
develop and numerically solve a modified version of the 
previously presented simulation of a simple physical model 
for this quicksand/collapsing soil. Here, we focus on the 
penetration behavior of an intruder. 



Experiments and Model 

Our physical model is inspired by in situ measurements 
performed with a specific type of natural quicksand at the 
shore of drying lagoons located in a natural reserve called 
Lengois Maranhenses in the North-East of Brazil [2011211 
22 ( 23] . Cyanobacteria form an impermeable crust, giving 
the impression of a stable ground. After breaking the crust 
a person rapidly sinks to the bottom of the field. We mea- 
sured the shear strength of the material before and after 
perturbation and found a drastic difference. Our measure- 
ments indicate that the quicksand is essentially a collaps- 
ing suspension with depth independent shear strength. Af- 
ter the collapse, it becomes a soil dominated by the Mohr- 
Coulomb friction criterion for its shear strength. The ma- 
terial undergoes a cross-over from a yield stress material, 
i.e. a more fluid-like behavior to a Coulomb material, i.e. 
more solid-like behavior after the collapse. We would like 
to point out that the collapse of the metastable structure is 
irreversible, as opposed to quicksand described in Ref. |17| . 
In summary, the "living quicksand" studied here can be 
described as a suspension of a tenuous granular network 
of cohesive particles. If perturbed, this unusual suspen- 
sion can drastically collapse, promoting a rapid segrega- 
tion from water, to irreversibly bury an intruding object. 

To model the complex behavior found in the exper- 
iments we use a variant of contact dynamics, originally 
developed to model compact and dry systems with lasting 
contacts [24ll25j . The absence of cohesion between parti- 
cles can only be justified in dry systems on scales where 
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the cohesive force is weak compared to the gravitational 
force on the particle, i.e. for dry sand and coarser mate- 
rials, which can lead to densities close to that of random 
dense packings. However, an attractive force, e.g. due to 
capillary bridges or van der Waals forces, plays an impor- 
tant role in the stabilization of large voids [55] > leading to 
highly porous systems as e.g. in fine cohesive powders, in 
particular when going to very small grain diameters. Also 
for contact dynamics a few simple models for cohesive par- 
ticles have been established [27ll28l[26] . Here we consider 
the bonding between two particles in terms of a cohesion 
model with a constant attractive force F c acting within a 
finite range d c , so that for the opening of a contact a finite 
energy barrier F c d c must be overcome. In addition, we im- 
plement rolling friction between two particles in contact, 
so that large pores can be stable [56ll5W3"0ll31l[3"5] . 

In the case of collapsing "living quicksand" our cohe- 
sion model has to be modified. One also has to take into 
account the time necessary for bonds to appear, i.e. dur- 
ing relatively fast processes new bonds will not be formed, 
whereas during longer times bonds are allowed to form 
at a particle contact. Finally, gravity also cannot be ne- 
glected in the model since the particle diameter is usually 
well above the micron-size. For simplicity, however, the 
surrounding pore water is not explicitly considered but 
only taken into account as a buoyant medium, reducing 
the effective gravity acting onto the grains. Disregarding 
the interstitial fluid motion keeps the model as simple as 
possible and nevertheless able to reproduce the main ex- 
perimental observations [21]. The details of the collapse, 
however, may be influenced by the flow field of the sur- 
rounding fluid |191l3"5] . Testing the influence of the fluid 
by introducing a viscous drag onto the grains considering 
water and the typical grain sizes from the experiments 
showed no significant difference [2"2"ll2"3] . As previously 
discussed, we justify the cohesive bonds in this case as 
being mediated by the bacteria living in the suspension. 

Summarizing, we use the following contact laws for the 
simulations presented in this paper. Perfect volume exclu- 
sion (Signorini condition) is assumed in normal direction, 
where a cohesive force is added as described above. In 
tangential direction Coulomb friction is applied (Coulomb 
law). Additionally, for the contact torques rolling friction 
is applied as a threshold law similar to the Coulomb law. A 
detailed description of this model can be found in Ref. [26] . 
For the simulation of the collapsing living quicksand the 
formation and breakage of cohesive bonds has to be con- 
sidered additionally. Our physical model is validated with 
the real data obtained from the situ measurements, specif- 
ically by comparing the shear strength behavior showing 
a drastic change when soil was perturbed and the pene- 
tration behavior described in more detail in Refs. [2TJll2"Tl 
[22]. 
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Results 

The penetration behavior of our collapsing "living quick- 
sand" has been studied previously for one specific density 
[20,21,22,23 and is here only briefly summarized. The 
penetration causes the partial destruction of the porous 



network and the subsequent compaction of the disassem- 
bled material. We observe the creation of a channel which 
finally collapses over the descending intruder. At the end 
of the penetration process the intruder is finally buried un- 
der the loose debris of small particles. Furthermore, our 
simulations indicate that in the worst condition, when the 
cohesive force is completely restored, one could need a 
force up to three times one's weight to get out of such 
morass |22j . In this paper we will focus on the penetration 
process for varying densities while keeping the cohesive 
force constant. 

When creating fragile structures by ballistic deposition 
and settling due to gravity the density is determined by 
the strength of the cohesive force. Here we will show a way 
of generating fragile structures with determined density 
(within a given range). The particles are deposited ballis- 
tically, and settle due to gravity. After the settlement of 
all particles, the cohesive forces between them are tuned 
to the point in which a barely stable structure of grains is 
assured. This results in a tenuous network of grains, like in 
a house of cards (see fig. [T^,) , giving the maximal possible 




Fig. 1. Maximal density of a configuration after settling with- 
out change (a) and minimal density after removing all possible 
loose ends (b). 

density for this specific configuration. In the case shown 
the volume fraction is 0.441. For different configurations 
(different seed for the random number generator), used for 
later averaging, this maximal volume fraction has values 
ranging from 0.428 to 0.451. 
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For averaging we need configurations with the same 
initial density. How can we create such structures? Due 
to the ballistic deposition there are many loose ends in 
the structure, which do not carry any load. When elim- 
ination those loose ends successively one can reduce the 
density to a given value. This eliminating procedure will 
be briefly discussed in more detail. Particles with only one 
contact are not contributing to the force network carry- 
ing the weight of the material. Thus, these particles can 
be eliminated without leading to a collapse of the struc- 
ture. We check for each particle (in random order) if they 
have only one contact in which case we eliminate the par- 
ticle. After we went through all the particles we start 
again checking all particles, as there will be new parti- 
cles with coordination number one due to the elimination 
of a former neighboring particle. This can be done un- 
til one reaches a state where no loose ends are present 
any more. This defines the minimal possible density for 
a specific configuration as shown in fig. [lb, with volume 
fraction 0.344 for this specific case. For different config- 
urations this minimal volume fraction ranges from 0.327 
to 0.35. This process can be interrupted when a desired 
density is reached. 
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Fig. 2. Penetration depth here measured by the weight w 
of all grains above the intruder, i.e. larger vertical position, 
normalized by the total weight wt) depending on the force Fd 
acting on the intruder (normalized by cohesive force F c ) for 
different initial densities of the material. (Connecting lines are 
shown here for different curves to be better distinguishable.) 
Inset: Scaling the force axis by the volume fraction of the initial 
configuration leads to a relatively good collapse for low values 
of Fd- For higher forces the curves show small deviations due 
to inertial effects. 



Fig. [5] shows the penetration depth depending on the 
applied force for different densities of the initial config- 
uration. Reducing the density the threshold, needed for 
pushing in the intruder, reduces, as one expects. Scaling 
the force with the density, i.e. dividing the z-axis by the 
volume fraction leads to a data collapse for small forces 
only (Fig. [5J inset). In this regime where the intruder is 
relatively slow, the density fully determines how deep the 
intruder can be pushed in. For higher forces the collapse 
gets worse due to inertial effects. In this region a better 
collapse could be achieved when scaling with the square 
of the density. The influence of inertia can be understood 



as follows: Structures with larger densities strongly hinder 
the motion of the intruder. Thus, inertial effects are less 
effective for these structures, leading to lower penetration 
depth within these structures. 

When looking in more detail at the dynamics of the 
penetration process it can be seen that, after an initial 
phase, the velocity fluctuates around a constant value (not 
shown here) before finally decelerating to zero at the final 
intruder position, i.e. at the bottom, or at intermediate 
values. For different forces applied on the intruder these 
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Fig. 3. Relation between applied force and measured average 
velocity when pushing downwards at one specific density (y = 
0.432). A parabolic fit represents the data quite well. 



average velocities can be measured. The relation between 
the applied force and the measured velocity can be fitted 
by a quadratic function (fig. [3]), where the minimum of 
the parabola is very close to u avg = 0. Alternatively, this 
minimum can be fixed at w avg = leading to an almost 
identical fitting curve. Summarizing, we found: 



Fa - F thr oc v a 



(1) 



In this case the value for F t hr/F c is around 0.42 which is 
in accordance to the value one could estimate from the 
force dependence of the relative weight w/wt above the 
intruder after penetration. Let us consider the applied 
force to equalize the drag force on the intruder exerted by 
the surrounding "complex fluid" . The results indicate a 
yield-stress fluid behavior as also found in shear strength 
measurements f2"2"l[2"B"] . Future rheological measurements 
could define the specific rheological model. A square ve- 
locity dependence of the drag force usually implies that 
inertia effects are important. Here, the intruder has to ac- 
celerate the grains as they arc pushed downwards within 
the compaction process. On the other hand, the force is 
also needed to break cohesive bonds. The fact that we do 
not find a viscous like behavior (linear in velocity) agrees 
with the non-reversibility of the system. 

A link can be drawn to the pinning-depinning transi- 
tion for a force driven interface where the driving force Fd 
has to overcome a threshold Fthr leading to a power law: 



v oc (F d - F t 



thrj 



(2) 



The former quadratic behavior previously observed would 
lead to an exponent 9 = 1/2. The results obtained for 
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Fig. 4. When plotting the average velocity versus the ap- 
plied driving force Fd minus a threshold force Fthr, a power 
law fit has an exponent of around 1/2 for the three different 
investigated volume fractions. 

different densities are shown in fig|4j Note that one has 
to determine -Fthr by adjusting to the optimal power law 
behavior, then a power law fit estimates the exponent. 
This procedure gives some additional information about 
the error of F trir and 8. In the cases presented here varying 
-fthr by about 5% still showed relatively good power laws 
leading to changing the exponent by less than 5% which 
also serves here to determine the error of the exponents 
(instead of the lower value obtained by the statistical anal- 
ysis by the regression). The values for the exponents are 
0.5 ± 0.02 (for v = 0.43 and v = 0.35) and 0.49 ± 0.02 (for 
v = 0.4). Obviously for all densities the exponent agrees 
very well with 1/2, suggesting a quadratic dependence on 
the drag force acting on the intruder. 
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Conclusion 

We investigated the density dependence of the penetra- 
tion behavior of a model for collapsing "living quicksand" . 
We could achieve a data collapse for small forces when 
plotting the penetration depth depending on the applied 
force divided by the density. For higher forces or penetra- 
tion depth incrtial effects are important and the scaling 
is less pronounced. During the penetration process the 
intruder velocity fluctuates around a constant value ex- 
cept for the very initial acceleration and the final decel- 
eration. This constant velocity shows a power law with 
exponent 6 = 1/2 as function of the driving force minus 
a threshold force. This means that the drag force on the 
intruder shows a quadratic velocity dependence. Here the 
drag force is the driving force shifted by a threshold force. 
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